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We consider a soluble model of multifragmentation which is similar in spirit to many models 
which have been used to fit intermediate energy heavy ion collision data. We draw a p-V diagram 
for the model and compare with a p-V diagram obtained from a mean-field theory. We investigate 
the question of chemical instability in the multifragmentation model. Phase transitions in the model 
are discussed. 
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' Statistical models of multifragmentation have long been used to explain data from heavy ion coUisons. Such a 
, model was first invoked for Bevalac results [1,2] and similar physical ideas but with many substantial variations 
'"5 ■ were subsequently used for intermediate energy heavy ion collisions [3-5]. In this work we consider a model of 
multifragmentation, variations of which have found many applications in the literature [6-11,13]. Thermodynamic 
properties of a simpler version of this model have also been discussed [15-17]. The model we use here has two 
kinds of particles but no Coulomb interaction. Throughout the rest of this paper we will refer to this model as the 
thermodynamic model. Typically the number of particles in our model is 200 although we also use systems containing 
as many as 1000 particles. While we could have easily included a Coulomb interaction term, our objective here is 
I different. The aim here is to test if because of two kinds of particles two features which have been discussed widely in 
^\ ' recent literature (from studies in mean field theory) persist in the thermodynamic model. These features are chemical 
instability (analogous to mechanical instability) and first order transition turning into second order. We therefore need 
CN to highlight some features which are present both in the thermodynamic model and in mean-field theories when mean 
field theories are applied to intermediate energy heavy ion physics. Typically mean-field theories use homogeneous 
infinite matter (hence no surface energy terms) and no Coulomb interaction. Finite systems with Coulomb and surface 
' terms have also been included [14] in mean-field models but this makes discussions more complicated and we want to 
'"^ , stay at the simplest level. As shown in [6] surface energy terms play an important role in a thermodynamic model and 
^ ■ are included. Moreover, since we will concentrate on two component systems, symmetry energy terms are included 
^ ', as they are also in mean field theories. 



II. THE THERMODYNAMIC MODEL 

The thermodynamic model has been described in many places [6,7,11]. For completeness and to enumerate the 
parameters we provide some details. 

Assume that the system which breaks up after two ions hit each other can be described as a hot equilibrated nuclear 
system characterized by a temperature T and a freeze-out volume V within which there are A nucleons {A = Z + N) . 
The partition function of the system is given by 



Q-^--En^ (2-1) 
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Here ni_j is the number of composites with proton number i and neutron number j and LUij is the partition function of 
a single composite with proton, neutron numbers i,j respectively. The sum is over all partitions of Z, N into clusters 
and nucleons subject to two constraints: J2i j ^^ij — ^ ^^'^ Si j — ^- These constraints would appear to make 
the computation of Qz,n prohibitively difficult but a recursion relation exists which allows the computation of Qz,n 
quite easy on the computer even for large Z ov N [12]. Three equivalent recursion relations exist, any one of which 
could be used. For example, one such relation is 

Qz,7i ^ ^ '^^i,jQ z — i,n—j (^■^) 
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The average number of particles of the species i,j is given by 



Qz-i,N-j 
Qz,N 



(2.3) 



All nuclear properties are contained in ujij. It is given by 



^(24i+j)mT)3/2 
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Here Vf is the free volume within which the particles move; Vf is related to V through Vf = V — Vex where Vex is the 
excluded volume due to finite sizes of composites. This is the only interaction between clusters we try to simulate. 
Thus the thermodynamic model is not an exact description of the system considered here but another approximation 
to it which has some interesting features that we hope to show. This restricts the validity of the model to low density 
(i.e., large V). Further, we take Vex to be fixed, independent of multiplicity. In reality. Vex should depend upon 
multiplicity [18]. We take it to be constant and equal to Vq = A/ po where po is the normal nuclear density and A is 
the number of nucleons of the disassembling system. As in previous applications, we restrict the model to freeze-out 
densities less than p/po = 0.5 that is F > 2Vo. The factor qij^int is the internal partition function of the composite. 
Define a = i + j. Then 



Here W = 15.8 MeV, a = 18 MeV, s=23.5 MeV and e = 16.0 MeV. The reader will recognise the volume, surface 
and symmetry energy of the cluster i,j and the contribution to the internal partition function from excited states 
in the fermi-gas formulation. For a(= i + j) > 5 we use this formula. For lower masses we simulate no Coulomb 
case by setting the binding energy of ■'He=binding energy of "^H and binding energy of ^Li=binding energy of ^H. In 
the weight of eq 2.4 we have not included a Fisher droplet model r which is a power law prefactor that is important 
around the critical point. Away from a critical point, exponential terms dominant the weight and this is the region 
we study in this paper. Such a term can be included, but our main focus is on the role of the symmetry energy in 
two component systems and related questions of chemical instability. 

For a given a, what are the limits on i(or j = a — i)l This is a non-trivial question. In the results we will show, 
we have taken limits by calculating the drip lines of protons and neutrons as given by the above binding energy 
formula. Limiting oneself within the drip lines is a well-defined prescription, but is likely to be an underestimation 
since resonances show up in particle-particle correlation experiments. On the other hand, for a given a, taking the 
limits of i from to a is definitely an overestimation. 

There is another consideration which restricts the validity of the model. We have assumed (Eq.(2.1)) that the 
standard correction rii^jl takes care of antisymmetry or symmetry of the particles. In the range T > 3 MeV and 
p/po < 0.5 this is usually true. At low temperatures where one might apprehend the usual correction to fail, it survives 
because many composites appear, thus there is not enough of any particular species to make (anti)symmetrisation 
an important issue. At much higher temperature the number of protons and neutrons increase but as is well-known, 
the n! correction takes the approximate partition function towards the proper one at high temperature. We define 
y = Z/{Z + N) where Z and N are the total proton and neutron numbers of the disintegrating system and the 
theory works even at low temperatures if y is in the vicinity of 0.5. But for example, at T=5.0 MeV and y=0 (neutron 
matter) , this is a terrible model. Now the number of neutrons is large and the temperature is not high and Fermi-Dirac 
statistics must be enforced. This was studied quantitatively in [19]. In our applications of the thermodynamic model 
we will confine y to be between 0.3 to 0.7 and T >3 MeV. This is indeed not very restrictive since this encompasses 
the drip lines and so the model, which was devised for intermediate energy heavy-ion collisions, will be applicable. 
For ^^^Sn-|-^^'*Sn coUisions, a much studied case, the value of y is 0.4. 

By equation of state (EOS) we mean p—V diagrams for fixed temperatures. This can be obtained by exploiting the 
equation p = T^^^#^. From eqs (2.1) and (2.3), this reduces to p = yj{^{nij) — 1) where -1 within the parenthesis 
corrects for centre of mass motion. The simplicity of this formula suggests that we just have a non-interacting gas 
with many species with V being replaced by Vf. But this is deceptive. In fact (which is the multiplicity) is 

not fixed but varies as a function of both temperature T and volume V thus this is not anything as simple as a mixture 
of non-interacting species. It is indeed interactions which make or break clusters to produce the final equilibrium or 
statistical distribution of fragments. In this sense interactions are also included. 

Since ours is a canonical model, we do not need the chemical potentials /ip (proton chemical potential) and /i„ but 
we compute them anyway from the relation fi = (§^)y,T- We know the values of Qz,n,Qz-i,n and Qz,n-i- Since 
F is just —TlnQ, we compute fip from /Xp = —T{lnQz,N — lnQz-i,N) and /i„ from —T{lnQz,N — lnQz,N-i)- Indeed, 



qi,j,int = exp [Wa - era- 
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the grand canonical version of the thermodynamic model we are solving has been known for a long time in heavy ion 
collision physics [2] . There the /ip and /i„ arise naturally. We have checked that the grand canonical values of /ip and 
Hn are indeed very close to the ones we derive by exploiting the canonical partition functions whose values we know 
numerically. Throughout this work, whenever we plot ji's we have obtained the values from a canonical calculation. 
One might think that since our model has many species there should be many /x's but in fact all /i's can be expressed 
in terms of only fip and /i„. Since our model is based solely on phase space, chemical equilibrium is in fact implied. 



III. A MEAN-FIELD MODEL 



We want to contrast the model above with mean- field theories. Our mean- field calculation uses the simplest 
model consistent with nuclear matter binding energy, saturation density, compressibility and symmetry energy for 
asymmetric matter. The potential energy density is taken to be 

V{Pn,Pp) = ^PnP, + ^{pI + pI) + -^^ (3.1) 
Po ^PO tJ + i Po 

Here po=-l6 fm"*^ and p„ and pp are neutron and proton densities and p = pn+pp- The dimensionless constant a and 
Au,Ai, B (all in MeV) are chosen to reproduce nuclear matter binding at 16 MeV per particle, saturation density at 
compressibility at 201 MeV and symmetry energy at 23.5 MeV. The energy per particle (including kinetic 

energy) at T = is 

El A = a/-^ + ^{pl + pD + -^{^r + 22.135 x [^(^)2/3 + ^(^)2/3] (33) 
PPa 2ppo " a + 1 Po P Po) P Po) 

The values of the constants are: a = 7/6; A„ = -379.2 MeV; Ai = -334.4 MeV; B =303.9 MeV. 
The Hartree-Fock energy of an orbital is given by 

e = pV2m + Auipu/Po) + A,{p,/po) + S(p/po)" (3.3) 
The value of pp is found by solving for a given pp and (3 = l/T 

h^Jo exp(/3(ep-/.p)) + l ^ • ' 

Similarly /i„ is extracted from p„. The pressure has contributions from kinetic energy and the potential energy. The 
contribution from the kinetic energy is calculated from well-known Fermi-gas model formula. Contribution to pressure 
from interaction is 

Pskyrme = PnPp + -^[{Pn) + (Pp) J + — T-^i — )^ (3-5) 

po 2po cr + i- po 



IV. EOS IN THE TWO MODELS 



Fig. 1 compares the p — p diagrams at constant temperature for the two models. We restrict the value of y = 
Pp/{Pp + Pn) between 0.3 and 0.5 and p/po between and 0.5 because outside these ranges the validity of the 
thermodynamic model is significantly reduced although the mean-field model has no obvious limitations. In Fig.l two 
curves for the mean-field model arc shown for each temperature and y; one is the straightforward p ~ p diagram and 
the other is with Maxwell construction which eliminates the spinodal instability region. While the straightforward 
p — p diagram in the mean field theory is widely different from the the p — p diagram of thermodynamic model, the 
p — p diagram with Maxwell construction is much closer specially at temperature 7 MeV. The p — p diagrams at T=10 
MeV in the two models are not that close. At least part of the reason is that the thermodynamic model is drawn for 
exactly 200 nucleons but the mean-field theory uses a grand canonical ensemble and hence is applicable to infinite 
systems only. We have checked that for a system of 500 nucleons the thermodynamic p — p diagram is closer to the 
Maxwell constructed mean-field p — p diagrams. 
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Fig.l: The EOS in the two models at different tem,peratures 7 and 10 MeV (left and right panel respectively) and 
with different proton fractions. The dashed and dotted lines represent the EOS in the meanfield and the 
thermodynamic model and the solid line is the Maxwell construction in the meanfield model. 

We want to point out that regions of negative compressibility (dp /dp < 0) which are common in mean-field theory 

(Maxwell construction eliminates these) arc almost absent in the thermodynamic model (they arc present when plotted 
in an expanded scale, see Fig. 3) and one would be tempted to conclude that the thermodynamic model is a good 
lowest order approximation. The thermodynamic model includes all inhomogeneous distributions of matter from 
single nuclcons and light clusters with gas- like behavior to very large liquid- like clusters. This feature approximates 
the Maxwell construction incorporated into a mean field theory which splits the system into two parts with liquid and 
gas densities. The very small region of negative compressibilty left over has its origin probably in the finite particle 
number effect and not is an inherent error in the model. This is dealt with again in section VIII. Mean field theory 
descriptions of two component systems introduce new features into the description not present in one component 
systems. Specifically, a new variable has to be introduced, such as the proton fraction y, which can be different 
in the two phases. The liquid phase can have one value of y, with the gas phase having another value, while still 
maintaining the total number of protons and neutrons. With a new variable, the coexistence curve and instability 
curve of one component systems become surfaces in p.T, and y for two components. In mean field descriptions, 
the first order phase transition of one component systems become a second order phase transition in two component 
systems. Moreover, mechanical instability and chemical instability no longer coincide. We now turn our attentions to 
features associated with chemical instability in our model. 



V. ISOSPIN FRACTIONATION 



Isospin fractionation is a well-established experimental phenomenon [20]. If the disintegrating system has a given 
N/Z > 1 then, after collision, the measured rin/np ratio (where n„, Up are measured single neutron and proton yields 
respectively) is higher than N/Z. Similarly, the ratio of measured < ni^2 > / < "•2,1 > is higher than what one might 
expect from the N/Z ratio of the disintegrating system. This then implies that if there is a large chunk left after the 
breakup it must have a. n/z ratio lower than original N/Z since the total number of neutrons and protons must be 
conserved. If we characterise n/z ratio etc. in terms of the parameter y we have been using, then if y source is less 
than 0.5, then y of the large chunk is greater than ysource and < Up > /(< rip > + < n„ >) is less than ysource- 

A priori, it would seem difficult to get this aspect out of a mean-field model but in a seminal piece of work MuUer and 
Serot have demonstrated how this might come about [21,22]. In mean-field theory, analogous to mechanical instability 
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(ip < 0) there appears regions of chemical instability, i.e., < (or > 0) when two kinds of particles are 
involved. One can avoid this unphysical region of chemical instability but then needs to consider splitting the system 
into two parts, each homogeneous but distinct from each other, one belonging to the liquid phase with higher y value 
and the other to the gas phase with lower y value. One consequence of this is that the phase transition takes place 
neither at constant pressure nor at constant volume and what would have been a first order phase transition, becomes 
a second order phase transition. 
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Fig. 2: Example of is o spin fractionation in the thermodynamic model, y of the largest cluster (top panel) is plotted in 
the top panel. This y is larger than the y value of the whole system. The lower part of the figure shows that the gas of 
nucleons is very rich. For example for y=0.4 and T=7-10 MeV, {un) w lO(np) while for T=4 MeV, {un) « lOO(np). 

In the thermodynamic model, isospin fractionation happens naturally. In general, the model has, as final products, 
all allowed composites, a,b,c,d.... where the composite labelled a has ya = ia/iia + ja) where ia,3a is the number 
of protons and number of neutrons respectively in the composite a. The only law of conservation ]s Z = ia x Ua 
and N = '^^ja x JT-a- So a large chmik can exist with higher y than that of the whole system and populations of 
other species can adjust to obey overall conservation laws. Whatever partition lowers the free energy will happen. 
The thermodynamic model is dramatically different from mean field models. The most significant difference is that in 
the thermodynamic model, if we prescribe that dissociation takes place at p/ pa=0.3 we still have only clusters with 
normal nuclear density and properties and also nucleons. It is just that there are empty spaces between different 
clusters and nucleons in the region of dissociation. But in mean field models p/po=0.3 will imply that the nuclear 
matter is uniformly stretched to this density. While this can happen as a transient phenomenon such as in transport 
calculation, whether this can also exist as an equilibrium situation is highly questionable. 

An example of isospin fractionation in the thermodynamic model is shown in Fig. 2. The formalism developed in 
section II can also be extended to calculate the average number of nucleons and protons (or neutrons) in the largest 
cluster. For brevity we do not write down the formulae here but these are straightforward extensions of Eqs.(2.7) and 
(2.8) given in [6]. Fig. 2 shows results from such a calculation. If y of the disintegrating system is less than 0.5, the 
y value of the largest cluster is larger than that of the source. Correspondingly, {{rip}) / {{np) + (n„)) is much smaller 
that y source- [It should be mentioned that the number of protons and neutrons will be augmented from decays of hot 
composites, so what is plotted in Fig. 2 is not what will actually be observed in experiments]. Further the isospin 
fractionation happens whether the dissociation takes place at constant volume or constant pressure. 

In this and many other aspects, the thermodynamic model is very similar to the Lattice Gas Model (LGM) with 
isospin dependence. For an accurate solution of LGM one has to give up the mean-field approach and obtain results 
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by Montc-Carlo simulation. Here also many composites are produced with many different y values [23,24]. Isospin 
fractionation happens naturally [24]. 



VI. INSTABILITY IN THE THERMODYNAMIC MODEL 

Fig. 1 shows that compared to the mean-field model, regions of mechanical instability with < nearly disappear 
in the thermodynamic model. In an expanded scale, they are more readily seen (Fig. 3) where we have drawn p — p 
diagram for a constant temperature T=7.0 MeV but different y's. 




Fig. 3: EOS (p — p) in the thermodynamic model at T=7.0 MeV but for different y 's. 

We clearly have some regions of mechanical instability. Chemical instability implies (^^)p,T < 0. We investigate 
that now. At T=7 MeV, we have drawn Hp (and /i„) at four pressures (Fig. 4). To get an understanding of the 
behaviour, we need to also look at Fig. 3. At the lowest pressure shown, p = 0.02 MeV fm~^ (Fig. 4), the horizontal 
constant pressure curve cuts the isothermals (Fig. 3) at the low density side only (between A and B) and pp rises 
monotonically between y = 0.3 and y = 0.5. The next constant pressure curve, at p = 0.025 MeV fm~^ (Fig. 4) cuts 
all isothermals (Fig. 3) at low density side {p/po < 0.1) between C and D and a few isothermals at higher density 
side. Between C and D, y increases as does pp. The points marked D and E have the same values of p and T but 
very slightly differing values of pp. As we move to the right from E along the line p=0.025 MeV fm,~^ the value of 
y drops as also the value of Pp. We forego describing graphs at other pressures but the figure shows there is a very 
small region where is negative ( Fig. 4, p = 0.035 MeV fm~^). The not so obvious feature is the appearence of 
two branches in both pp and Pn (i.e., for example, the p = 0.025 MeV fm~^ curve). The two branches would merge 
for a Van der Waals fluid and thus the appearence of two branches indicates departure from a Van der Waals fluid. 



VII. COMPARISON WITH VAN DER WAALS FLUID 



For a Van der Waals fluid with Maxwell construction, the following behaviour will be seen [25] as we move along 
an isothermal ina p — p plane, provided we are below the critical temperature. If we start with very small density we 
are in the gas phase. As the density rises, the chemical potential changes till it reaches the coexistence region. In this 
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region /igas = I^Hq and as the density changes the chemical potential remains unchanged but more particles change 
from the gas phase to the liquid phase. This remains the situation till a high density is reached when all the particles 
are in the liquid phase. 





Fig. 5: fj, = yup + (1 — y)A*ra as a function of p/ po for different y 's(left panel). The temperature is 7 MeV. In the right 
panel the behaviour of p, is shown as a function of pressure p for y = 0.5 and T = 7 MeV. 

The situation in the thermodynamic model is depicted in Fig. 5. The model is inapplicable to high density situation, 
so we cut-off p at p/po = 0.5. We have two chemical potentials fXp and pn and to compare to a Van der Waals fluid, 
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we consider the combination fi = yfip + (1 — y)j-J,„ (for y=0.5, /ip = /i„ any way). At T=7.0 MeV we show in the left 
panel of Fig. 5 the behaviour of /i at j/ = 0.5, 0.4 and 0.3. For the y=0.5 curve we also depict schematically what the 
behaviour would have been if we had a Maxwell-corrected Van der Waals fluid. From some point A, the chemical 
potential would remain unchanged (shown by a horizontal line ending at B which is the end-point of our density). A 
more familiar plot is n against pressure p for a fixed temperature. This is shown in the right panel for our model. For 
a Van der Waals fluid, the segment from AtoB would simply collapse to the point A. 

VIII. SPECIFIC HEATS IN THE MODEL 

In [6] where the thermodynamic model was first studied for phase transitions, it was pointed out that for a given 
density p, the specific heat per particle Cv/A tends to oo at a particular temperature when the particle number A 

tends to oo. Since Cy = (f^)y = ^(f^)v = ~'^{§^)v, a singularity in Cy signifies a break in the first derivative 
of F, the free energy and a first order phase transition. The model in [6] considered one kind of particle although 
binding energy, surface energy etc. were chosen to mimic the nuclear case. Wo see similar effect here when we take 
into account two kinds of particles explicitly (Fig. 6, see also [7]). The calculated Cy/A becomes progressively sharply 
peaked as A increases for all y values between 0.3 and 0.5. This behaviour of the specific heat is very different from 
that of mean field model of nuclear matter where the specific heat at constant volume varies smoothly from a low 
temperature Fermi gas to an ideal gas as T increases. In a thermodynamic model with fragmentation, this behaviour 
is modified by the surface energies that arise in the multifragmentation of the original nucleus into clusters of different 
sizes. The peak in the specific heat occurs at the point where the largest cluster suddenly disappears. This behaviour 
is nuclear boiling. 

Specific heat per particle Cp/A in the model has not been considered before. We notice there are regions with 
^ < (even though these regions are much less visible than in the mean field model) . Their presence might indicate 
finite particle number effects or it may be a shortcoming of the model. Calculations which are on going suggest it is 
particle number effect rather than an inherent problem in the model. These negative regions of compressibility can 
lead to negative values of Cp. This is a contentious issue at the moment and we intend to deal with these issues fully 
in a future publication. 
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Fig. 6: Cv/A as a function of temperature for systems of 200, 500 and 1000 particles with different proton fractions. 
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IX. SUMMARY 



We looked at several features of a thermodynamic model (which has seen many applications in data fitting) and a 

mean field model. Equation of state in mean-field theory has large regions of mechanical instability even for infinite 
systems and one needs to do a Maxwell construction to eliminate these. By contrast, the thermodynamic model 
has directly an EOS which becomes very flat with density and volume and this behaviour resembles a real system 
undergoing a first order phase transition. In mean-field models, when the system enters the region of instability, it 
fragments into pieces. This fragmentation is directly included in the thermodynamic model and this is the reason 
for relative flatness in the EOS. The cluster distribution readjusts itself with changes in y or p to maintain a nearly 
constant pressure. Isospin fractionation seen in experiments can be also obtained in the mean field model but it 
requires a bifurcation in the isotopic space. It also requires that during dissociation neither pressure nor volume 
remain constant. By contrast, isospin fractionation occurs naturally in the thermodynamic model and can happen 
either at constant volume or at constant pressure. Large differences between these two models also appear in the 
calculation of Cy- The thermodynamic model has a strong peak in Cy whose origin are surface energy terms in 
the multifragmentation process which is lacking in a mean field model of homogeneous nuclear matter. This peak is 
associated with the phenomenon of nuclear boiling and the sudden disappearence of the largest cluster. 
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